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Non-local electrostatic interactions associated with the finite solvent size and ion polarizability are 
investigated within the mean-field linear response theory. To this end, we introduce a field theoretic 
model of a polar liquid composed of linear multipole solvent molecules and embedding polarizable 
ions modeled as Drude oscillators. Unlike previous dipolar Poisson-Boltzmann formulations treating 
the solvent molecules as point dipoles, our model is able to qualitatively reproduce the non-local 
dielectric response behavior of polar liquids observed in Molecular Dynamics simulations and Atomic 
Force Microscope experiments for water solvent at charged interfaces. The present theory explains 
the formation of the associated interfacial hydration layers in terms of a cooperative dipolar response 
mechanism driven by the reaction of the solvent molecules to their own polarization field. We also 
incorporate into the theory the relative multipole/dipole moments of water molecules obtained 
from quantum mechanical calculations, and show that the multipolar contributions to the dielectric 
permittivity are largely dominated by the dipolar one. We find that this stems from the mutual 
cancellation of the first two interfacial hydration layers of opposite net charge for multipolar liquids. 
Within the same non-local dielectric response theory, we show that the induced ion polarizability 
reverses the interfacial ion density trends predicted by the Poisson-Boltzmann theory, resulting in a 
surface affinity of colons and exclusion of counterions. The results indicate that the consideration of 
the discrete charge composition of solvent molecules and ions is the key step towards a microscopic 
understanding of non-local electrostatic effects in polar solvents. 

PACS numbers: 03.50.De,05.70.Np,87.16.D- 



I. INTRODUCTION 

The precise determination of electrostatic interactions 
in the vicinity of charged molecules in water solvent is 
one of the biggest challenges in colloidal sciences. From 
the performance of energy storage devices [lH3] and wa- 
ter purification membranes [H [S] to the solubility of salt 
ions and polyelectrolytes in water [7], a wide variety 
of electrostatically driven nanoscale processes depend on 
the electrostatic potential behavior close to charged sub- 
strates. Water that mediates these electrostatic interac- 
tions being a strongly polar liquid, the evaluation of the 
electrostatic potential requires in turn a proper insight 
into the dielectric response of solvent molecules to the 
charged sources. Furthermore, the ordering of solvent 
molecules at charged mica surfaces revealed by Atomic 
Force Microscope (AFM) experiments [H] and Molecular 
Dynamics (MD) simulations [31 HI [TD] indicates that the 
electrostatic interactions in real systems are non-local. 
Therefore, a consistent formulation of non-local electro- 
statics is needed to understand nanoscale phenomena. 

Our limited understanding of the dielectric response of 
water mainly stems from the lack of a microscopic the- 
ory able to map from the molecular details of the solvent 
to experimentally accessible macroscopic observables. In 
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particular, dielectric continuum theories that bypass the 
charge composition of solvent molecules are unable to 
account for their coordinated behavior in the presence 
of charged sources. This approximation resulting in a 
uniform dielectric response characterized by a constant 
dielectric permittivity s{r) — has severe drawbacks. 
For example, in bulk polar liquids, the local Born theory 
that cannot account for the dielectric void around ions 
is known to strongly overestimate ionic solvation ener- 
gies [TT]- 

Macroscopic theories of non-local electrostatics based 
on a phenomenological dielectric permittivity function 
e(r, r') and providing better agreement with experiments 
have been proposed over the last three decades. Among 
several efforts in this direction, one can mention the 
seminal works from A. A. Kornyshev et al. that dealt 
with non-local effects on the solvation of ions in bulk 
liquids [m [T3] and the charge storage ability of metal- 
lic capacitors |14) . Different formulations incorporating 
dipolar correlation effects in a coarse-grained way have 
been also developed in order to improve over the local 
Born theory in bulk [11 and confined solvents jl^. 

In inhomogeneous liquids, the Poisson-Boltzmann 
(PB) formalism based on the dielectric continuum ap- 
proximation fails as well to describe the interfacial dipo- 
lar ordering effect observed in experiments and simula- 
tions. The first dipolar Poisson-Boltzmann (DPB) ap- 
proach able to account for the electrostatics of solvent 
molecules was introduced in Ref. [16] . The excluded vol- 
ume of ions and solvent molecules was later incorporated 
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to this formalism at the mean-field (MF) level of approx- 
imation |17] , and the DPB formalism was reconsidered in 
bulk liquids at the one-loop order [T5]. Different gener- 
alized PB approaches based on dielectric continuum but 
accounting for the multipolar moments of ions were also 
developed in Refs. [19l|20]. We have recently extended 
the DPB formalism of Ref. [TB] beyond MF level in order 
to investigate in inhomogeneous electrolytes surface po- 
larization effects on the differential capacity of low dielec- 
tric materials [2]- This extended DPB (EDPB) approach 
allowed us to improve the agreement of the PB theory 
with experimental capacitance data of carbon-based ma- 
terials in a significant way. However, these generalized 
approaches that treat the solvent molecules as point- 
dipoles cannot consider their extended charge structure. 
As a result, they yield exclusively a local picture of elec- 
trostatic interactions in charged systems. At this stage, 
one should mention the innovative works of Refs. [2TH23] 
that focused on the electrostatics of charges with finite 
extension at different approximation levels, though the 
solvent was again considered in these models within the 
dielectric continuum approach. 

In order to overcome these limitations, a microscopic 
polar liquid model accounting for the discrete charge 
composition of solvent molecules is needed. In this ar- 
ticle, we introduce a microscopic description of non-local 
electrostatic interactions in polar liquids. We first derive 
in Sec. |lT] the field theoretic model of the polar liquid 
composed of multipolar solvent molecules of finite size, 
and embedding polarizable ions modeled as Drude os- 
cillators. Then, we obtain from the saddle point of the 
model Hamiltonian a non-local PB equation, which is 
considered in the rest of the article within the linear re- 
sponse regime of polar liquids symmetrically partitioned 
around weakly charged planar interfaces. Sec. |III| intro- 
duces a mapping from the microscopic polar liquid model 
onto the macroscopic relations of non-local electrostatics. 



Finally, we investigate in Sec. IV dipolar correlation ef- 
fects, multipolar contributions to the permittivity of the 
liquid, and the effect of induced polarizability on the in- 
terfacial solvent and salt partition. The limitations of 
the present theory, its possible extensions and potential 
applications are thoroughly discussed in the Conclusion 
part. 



II. FIELD THEORETIC MODEL AND MF 
EQUATIONS 

This part is devoted to the derivation of the field the- 
oretic model for a polarizable ion gas of different species, 
immersed in a solvent composed of polar molecules with 
a linear charge distribution. A schematic presentation of 
the solvent charge geometry is given in Fig. [ija). Each 
solvent molecule consists of a rigid rod where n,. elemen- 
tary charges with valency Qi are distributed, with the 
index I running over the elementary charges on the sol- 
vent molecule. Our consideration of linear multipolar sol- 
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FIG. 1: (Color online) (a) Geometry of polarizable ions (left) 
and multipolar solvent molecules (right), (b) Charged inter- 
face. 



vent molecules is motivated by the possibility to incor- 
porate into this geometry the relative dipole/multipole 
moments of water molecules obtained in quantum molec- 
ular calculations [23]. This complication will be treated 
in Sec. IIVDI We also note that each solvent molecule is 
overall neutral, that is Qi = 0. Furthermore, the dis- 
tance of the charge Qi from the first charge is a;, and 
fli = corresponds to the origin of the molecule. Fi- 
nally, — B. = a„^u is the vector pointing the end of 
the molecule, with a„^ = a the total molecular size and 
the unit vector u is parallel to the oriented molecule. 

The composition of polarizable ions of p species is also 
displayed in Fig. [ija). Each ion of species i consists of 
two elementary charges of valency and separated by 
the distance b. The electroneutrality condition implies 
J2iPibisi + Ci) = 0. Moreover, the ionic polarizability 
associated with the deformation of the electronic cloud 
by the surrounding fields is considered within the Drude 
oscillator model p5j . 
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where bp^ is the variance of these oscillations associated 
with the ions of species i. It will be shown in Sec. |IV E| 
that 6pj is proportional to the induced polarizability a. 
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The canonical partition function for the system of po- 
lar molecules and polarizable ions coupled with electro- 
static interactions read 
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where is the number of solvent molecules, Ni the num- 
ber of ions for the species i, and Xrd and Xxi stand for 
the thermal wavelength of solvent molecules and ions, 
respectively. We also introduced the shorthand nota- 
tion V ~ ({xfc}, {a^}, {yij}, {bj}) for the configurational 
vector space, where y^j and x^: stand respectively for 
the spatial coordinates of the first elementary charge 
of ions and solvent molecules. Finally, we note that 
= {dk,(pk) denote the solid angle characterizing the 
orientation of the fcth solvent molecule. 

The interaction energy is composed of an electrostatic 
part and a wall contribution, i/(v) = Hei[y) + i/u,(v), 
where the electrostatic part reads 
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[Pic + Psc + 0-]r «c(r, !•') [Pic + Psc + 0-]r' , 

(3) 

with (T(r) the fixed charge distribution (see Fig. [l]^b)), 
and the ionic and solvent charge densities corresponding 
to the charge compositions in Fig. [l]^a) are respectively 
defined as 

Psc{'c) = ^'^Qi5{Y-yik-&i)- (5) 

k=l 1=1 

Furthermore, the Coulomb potential in Eq. (|3| is defined 
as the inverse of the following operator. 



(6) 



where e(r) is a spatially varying dielectric permittiv- 
ity. We also note that the bulk Coulomb potential 
reads v^^{r) = ^s/r, with the Bjerrum length in the 
air £b{^) ~ ^ I [47r£(r)/cBT] ~ 55 nm, e the elementary 
charge, and T — 300 K the ambient temperature. The 
inverse of the bulk potential is given by the inverse of the 
following operator, 

-r\r,r') = -^A5(r~r'). (7) 

In Eq. ([t]), eo stands for the permittivity of the air. We 
will consider in this work exclusively the case of a uni- 
form background permittivity £(r) = Eqi which implies 
t'c(r, r') = v\(v — r') and ^s(r) = Ib- We also note 



that in Eq. ([2]) , we subtracted from the total Hamiltonian 
the self energy of ions Ei = (e| -f cf) Vc{Q)/2 + eiavdh) 
and polar molecules Eg = bc(0) — Vc{a)\ in the air 
medium. 

Finally, the part of the Hamiltonian corresponding to 
particle-wall interactions read 
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where we introduced the general wall potentials 
^iiyiji^j) and Ws(r, a) respectively for ions and sol- 
vent molecules. These wall potentials will be used to 
restrict the phase space accessible to the particles, and 
also as generatrice functions in order to derive particle 
number densities. We also note that all energies will be 
given in units of the thermal energy fc^T, the dielectric 
permittivities in units of the air permittivity Eq, and the 
surface charges in units of the elementary charge e. 

In order to simplify the theoretical analysis of the sys- 
tem, one can pass from the density to the field representa- 
tion by performing a Hubbard- Stratonovich transforma- 
tion. The grand canonical partition function defined as 
Zg - EAr.>onLiEjv.>oe'''^'e'^™'^=^c takes the form 
of a functional integral over this fiuctuating electrostatic 
potential, Zq ^ jVcj) e'"^'^\ with the Hamilto nian func- 
tional 
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where we rescaled the ionic and solvent fugacities as A.; 
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e^' /X^i and A^ = 

In order to derive local number densities, we first 
split the solvent and ion wall potentials into two parts, 
W^,(r,a) = Wsi{r) + Ws2(.r,a) and W^,(r,b) = Wa{r) + 
Wi2{r, b). Passing now from the complex to the real elec- 
trostatic potential with the transformation (j){r) — i(j){r), 
the mean-field level ion and solvent number densities re- 
spectively follow by taking the functional derivatives of 
Eq. ([9| with respect to Wii(r) and W^iir), 

db 



Pip{r) = Ai 
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-ei0(r) — Ci0(r+b) 



(10) 
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In terms of the same real electrostatic potential, the MF- 
level saddle point equation SH[(f)]/5(j){r) = takes the 
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form of a generalized PB equation, 



A0(r) + AttIb 



0, (12) 



where we introduced the ionic and solvent charge densi- 
ties 



Pici'r) = Pib 
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and used the MF relations between the charge densities 
and fugacities A^e^* = pn, and A^e^^ = psb- These re- 
lations follow from the bulk limits of Eqs. ( 10 ) and (111. 
We also note that in this work, the particle and charge 
partition functions for ions and solvent molecules will 



be related to the local densities in Eqs. (lOl-(ll) and 
Eqs. (13)-(14) according to fc(r) — p{r)/pb. 



The dependence of the solvent charge densities in 
Eq. ( 14 1 on the values of the electrostatic potential at 
different points around z makes Eq. (12 1 a non-local 



Poisson-Boltzmann (NLPB) equation that embodies the 
non-local dielectric response of the polar liquid at the 
molecular level of precision. We also note that for 
ions with vanishing polarizability {hp = 0) and solvent 
molecules of dipolar geometry (see Fig.[2]ja)), by expand- 
ing t he argument of the potential in the exponenti al of 
Eq. ([l4| at the order o\a?), the NLPB equation ([l2| 
tends to the DPB equation derived in Ref. |16) . 

We will investigate in this article the MF theory of non- 
local electrostatic interactions for polar liquids in contact 
with a charged planar interface located at z = 0, and 
corresponding to a surface charge distribution ct{v) = 
—asd{z) with as > 0. The negatively charged wall splits 
the space accessible to the liquid into two regions z < 
and z > 0, with equal dipolar and ionic bulk concentra- 
tions on each side (see Fig.[l](b)). The rotational restric- 
tion for dipoles will be considered exclusively at the end 
of the Sec. |IV C[ and we will assume that the interface 
at z = is penetrable in the rest of the article, that is 
Wsir,fl) = Wt{r) = 0. We also note that due to the 
translational symmetry within the {x, y) plan, the elec- 
trostatic potential depends only on the separation from 
the wall at z = 0, i.e. (j){r) = (j>{z). 

Moreover, we will consider exclusively the linear re- 
sponse regime corresponding to weak surface charges. By 
expanding Eq. ( 12 ) at the linear order in the electrostatic 
potential 0(z), and passing from the azimuthal angle 9 
to the projection of the dipole orientation on the z— axis 
with the transformation = a cos 9, one obtains the 



following non-local differential equation, 
A4}o{z) + Att£b(j{z) - eu,Kj(/)o(z) - Kl(j){z) 



(15) 
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X [0o(z + b,) + Mz - b,) - 20o(z)] 



0, 



where — AniBj^i PibQi /^w ^^'^ ionic screening pa- 
rameter, with qi = ei + Ci the total charge of each ion of 
species i, and nj. = AniBPsb J2i Qf ^^'^ screening parame- 
ter associated with solvent charges. We introduced above 
the notation aim = ai — am for the separation distance 
between the charges I and m on the solvent molecule (see 
Fig. [ij , and also the dielectric permittivity in the bulk 
solvent medium that will be defined in Sec. Ill (see 
Eq. (|22|)). The Bjerrum length in the water is indeed 
related to the one in the air medium as iw = is/^w 
Furthermore, the naught in 0o(-2^) means that deriving 
Eq. (15 1, we neglected the rotational penalty for the sol- 
vent molecules in the region z < a. We finaly note that 
all numerical results will be obtained for monovalent ions 
qi = 1 (until Sec. IV E on polarizable ions), and the model 



parameters of the dipolar solvent molecules will be cho- 
sen as Q = 1 and a = 1 A, unless otherwise stated. 



III. MAPPING TO THE MACROSCOPIC 
FORMULATION OF NON-LOCAL 
ELECTROSTATICS 

We introduce in this part a mapping from the micro- 
scopic model of Eq. ( [T5| onto the macroscopic formula- 
tion of non-local electrostatics. This mapping will allow 
us to relate the effective permittivity of the polar medium 
to the polarization charges in the liquid. By defining first 
the kernel operator 

G-i(z,z') = Z^t^^5{z-z') + ^5{z-z') (16) 
+Psby^QiQm [ ^^S{z'-z-a^) 

J -aim 
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X [S{z' -z-b^)+ 5(z' -z + b^)- 25{z' - z)] , 
and using the definition of the Green's function 



dz'G-\z, z')G{z', z") = S{z - z"). 



(17) 



the linear NLPB equation ( 15 ) that can be reexpressed 

as 



dz'G-\z,z')Mz') = ct(z') 



(18) 
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can be inverted, and the solution expressed in the form 

^+00 



dz'G{z,z')(Js{z'). 



Using now the relations (16 1 and (17 1, the Green's func- 



tion can be derived in ID Fourier space as 



G{z - z') = 44 



dk 



cos [k{z — z')] 



(20) 



where we introduced the dielectric permittivity function 
in Fourier space, 



k'^ 

Sir is 
fc2 



sin (fca;^) 



(21) 
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We note that the second-third and the forth terms on the 



r.h.s. of Eq. (21) correspond respectively to the charge 



structure factor of solvent molecules and polarizable ions 
in Fourier space. The bulk permittivity introduced in 
Eq. ( 15 ) is precisely defined as the infrared (IR) limit of 
Eq. (21), = e(fc — >■ 0), and it is given by 



= 1 ^ QiQmofjn - Stt^s ^ Pibbl^e^Ct. 

(22) 



Plugging now the Green's function (20 1 into Eq. (19), 



the electrostatic potential follows in the form 
2 Z""" .. cos(/cz) 



00 (^) 



dfc 



T^P-iqiJo K,f + k^e{k)/Su 



(23) 



where /ii = l/{2TTqi£was) stands for the ionic Gouy- 
Chapman length, and the net electrostatic field E{z) — 
4''q{z) reads 



k sin(fcz) 
TTAii^j Jo + k^^{k)/sy 



dk 



(24) 



We note that the potential in Eq. (23 



characterized by 
a diffuse permittivity function is simi ar in form to Eq. 
(3.16) of Ref. [H] where a phenomenological dielectric 
permittivity function was used in order to investigate 
non-local electrostatic effects on the charge storage abil- 
ity of metallic capacitors 



Using Eqs. (17) and (20), one can show that Eq. (16) 



can be recasted in the form of a non-local electrostatic 
kernel 



G"i(z,z') = 



1 



47r£f 



'dze{z ~ z')dz' + ewKf5(z ~ z')] , 

(25) 

where the non-local dielectric permittivity function is 
simply the inverse Fourier transform of Eq. (21 ), 



e(z - z') = 5{z - z') -f AirlBXiz - z'), 



(26) 



with susceptibility function 



(19) X{z) 



(27) 



- (z + a;„)^sign(z -f- aim)} 

^sr^ , ( \z\ ^ ,. f \z\ \ hi(z) 

+2^p,b6p,e,c,|^Erfc ' ' ' ^ ' 



where the distortion energy of polarizable ions hi{z) is 
given by Eq. ([I]). The equation (27) is the key result 
relating the microscopic polar liquid model of Eq. ^ to 
the macroscopic formulation of non-local electrostatics. 
A simpler form for this susceptibility function will be 
given for the case of simple dipolar liquids with point 
ions in Sec. |IV A| and for polarizable ions embedeed in 
the dipole liquid in Sec. |IVE| . 

We now note that in terms of the non-local dielectric 
displacement field 



D{z) 



dz'e{z~ z')E{z'), 



(28) 



the equation (18) can be rewritten as 

- d,D{z) = 47r£BO-(z) - e^nlMz)- (29) 

Neglecting the ionic screening term in the region k^z <C 1 
and integrating Eq. ( 29 ) , one gets for the induction field 



D{z) — 27r^B(TsSign(z). 



(30) 



We now define the polarization field P{z) through the 
usual non-local dielectric response equation 



P{z) 



dz'xiz- z')E{z'). 



(31) 



In the next section where we investigate the simplest case 
of ions without polarizability in a dipole liquid, it will be 
explicitly shown that the function x{z — z') brings the 
contribution from individual dipoles to the total dielec- 
tric correlation function, with a characteristic correlation 
length of the same order as the molecular size. Further- 
more, by using Eqs. (26) and (30), the polarization field 



introduced in Eq. (31 ) can be related to the displacement 



and total electrostatic fields as 



P(z) = 



Diz)-E{z) 
4tt£r 



(32) 



Reexpressed in the form E{z) = D{z) — AirtsPiz)^ this 
well-known macroscopic relation can be interpreted as 
follows. In a polar medium in contact with a fixed charge 
source (e.g. the surface charge located at z = 0), the local 
field experienced by a test ion at the position z is the 
superposition of the induction field D[z) generated by 
the surface charge (i.e. the field in the air medium), and 
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the reaction field induced by the polarizable molecules 
in response to this induction field. The reduction of the 
latter by the polarization field is the so-called dielectric 
screening effect. 

Furthermore, we note that in the same region HiZ <^ 1, 
the electrostatic field in Eq. ( 24 1 becomes 



E{z) 



QifiiEcsiz) ■ 



where we introduced the local effective dielectric permit- 
tivity function 



dk sin(fcz) 

T i{k) ■ 



(34) 



fnjecting now the relations (30 1 and (|33[) into Eq. (32), 



one obtains a relation between the macroscopic dielectric 
permittivity function and the polarization field. 



P(z) 



1l 
2 



sign(z) 



1 



(35) 



Moreover, by substituting the relations ( |3l[ ) and (32 1 into 
Eq. (29), and using the non-local PB equation (15), one 



obtains an expression relating the variations of the po- 
larization field to the polarization densities. 



dP{z) 
dz 



= Psbksc{z) + ^Pibkl,'}{z), 



(36) 



Thorough the expression (39), one can see the in- 
verse dielectric permittivity e^'-{z — z') as the dielec- 
tric correlation function containing the whole informa- 
tion on the polarizability of the solvent. Comparing the 
inverse permittivity in Eq. ( 40 ) with the effective permit- 
tivity Eq. ( [34| , we find that both functions are related 
e-i(z) = [2ecs{z)r\ With the use of Eqs. 



as 



(33) and (36 ), this relation shows that the inverse permittivity 



is related to the normalized charge densities associated 
with the polarizable molecules by the simple equation 



e-\z)=Siz)-^ 



Psbksc{z) + ^ Pifcfcpi (2 



(41) 



Injecting this relation into Eq. (39), we obtain a relation 



that expresses the modification of the induction field by 
the polar liquid as the convolution of the former with the 
charge density of polar molecules over the whole space. 



E{z) = D{z)^ 



dz' [psbksciz - z') (42) 



-Y,P^bk^Kz-z') D{z'). 



where the solvent charge partition function Eq. ( 14 ) takes 
in the linear potential approximation the form 



ksc{,z) 



Y^Q^Mz) 



(37) 



dflz 
2a/„ 



■(f)o{z + a^), 



and we also introduced the part of the ion charge parti- 
tion function associated with the ionic polarizability (i.e. 
the sixth term on the Ihs of Eq. (15)), 



+00 



^-h.{b,) 



(38) 



X [Mz + b,) + Mz - bz) - 20o(^)] • 
We now note that the inverse permittivity function that 



allows us to invert Eq. ( 28 ) as 



/oo 
dz'e-^{z- z')D{z'), 
-00 



is given by 



e-\z) ^ 



1 r°° dkcos{kz) 
nJo 1 + Msxik) ' 



(39) 



(40) 



where the Fourier transformed susceptibility is de- 
fined through the relations (21) and (26) as x(fc) = 
[e(fc)-l]/(4^£B). 



By considering now the explicit form of the induction 
field (30) in Eq. ( [42| ), one finally gets with the use of 
Eq. ( 33 ) a relation between the local value of the dielec- 



tric permittivity Scs{z) and the integrated polarization 
density. 



£c«{z) 



= 1 - 



Psbk, 



kp} (z ) 



(43) 

we used the reflection 
symmetry of the densities with respect to the interface. 



for z > 0. Deriving Eq. (43) 



i.e. ksci-z) = ksc{z) and kp}{-z) ^ k'^i{z). 

According to the relation ( [43| , in the linear response 
regime, local deviations of the dielectric permittivity 
from the permittivity of the air are related to the ac- 
cumulated polarization charge of solvent molecules and 
ions between the considered point in the liquid and the 
charged plan. As a result, the effective permittivity tends 
on the surface to the permittivity of the air, that is, the 
immediate vicinity of the interfacial area is character- 
ized by a dielectric void. The manifestation of this pe- 
culiarity absent in local electrostatic theories O [TBI Hi] 
but present in MD simulations ^3, 9, 10, and AMF ex- 
periments 18] indicates that the proper consideration of 
the extended charge structure of solvent molecules in our 
model is the key ingredient to recover the correct dielec- 
tric response of the water solvent. This point will be 
elaborated in further detail in the following parts. 
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FIG. 2: (Color online) Charge composition of (a) linear 
dipoles, (b) quadrupoles , (c) octupoles, and (d) polarizable 
ions considered in Sec. IIVEI 



IV. NON-LOCAL ELECTROSTATIC EFFECTS 
IN DIPOLAR LIQUIDS AT CHARGED 
INTERFACES 



We investigate in this part non-local electrostatic ef- 
fects on the dielectric response of the liquid, and the sol- 
vent and salt partitions at the charged interface. The 
geometry of the dipolar solvent molecules composed 
of two point-charges with valency ±Q is presented in 
Fig. [2]ja). For this specific solvent charge geometry and 
non-polarizable ions (i.e. bip = 0), we will first enlighten 
in Sec. |IV A| a collective dielectric response mechanism 
at the scale of single dipoles. Then in Sec. |IVB[ we will 
show that the same mechanism is responsible for the for- 
mation of successive hydration layers around the charged 
surface, which will be shown in Sec. |IVC] to explain the 
characteristic shape of the transverse permittivity pro- 
files observed in MD simulations |21 13 [10] ■ At the next 



step in Sec. IV D we will estimate the contribution of the 
multipolar moments of water molecules to the dielectric 
permittivity of the liquid. Finally, the effect of induced 
ion polarizability on the salt partition at charged inter- 
faces will be discussed in Sec. IIV El 



A. Collective dielectric response mechanism 

We examine in this part the non-local dielectric screen- 
ing mechanism discussed in Sec. |III| in terms of individ- 
ual dipole interactions. To this end, we note that for 
the dipolar charge composition in Fig. [2l^a), the dielec- 
tric susceptibility function defined in Eq. ( |27p is given by 
x{z — z') — 2pQPsb/aC{z — z'), with the dipole moment 



Pq = Qa and the adimensional susceptibility 

Ci{z -z') = Ul- ^^)^ ia-\z- z'\) . (44) 



A. The susceptibility function (44) shows that the polariz 



ability of the solvent liquid takes place over a finite region 
determined by the solvent molecular size. 

We now note that by expanding the denominator of 
the integrand of Eq. (40 ) in powers of K^a, and using the 
convolution theorem in Fourier space, the inverse permit- 
tivity function can be rewritten in the form of a geometric 
series, 

e-\z - z') = 5{z -z') + - V(-l)" («sa)'" Cr,{z - z'), 



n>l 



(45) 



where the high order correlation functions for n > 1 are 
related with the susceptibility in Eq. (44 1 by the recur- 
rence relation 



CJz- 



-C,{z^z")C,,_,{z" ~z'). (46) 



One notices that the first term on the rhs of Eq. ( 45 1 cor- 



responds to local dielectric correlations, and the second 
term of order O ((kso)^) that introduces the non-local di- 
electric response extends the range of these correlations 
by one molecular size a. It is also shown in Appendix \K\ 
that at the same order O ((ksc)^), the function C{z) is 
proportional to the charge density of a single dipole in- 
teracting with the charged surface in the air medium, i.e. 
ksc{z) — AQa/ p,sC{z), where we introduced the dipolar 
Gouy- Chapman length 



(47) 



Thus, the non-local contribution to the dielectric cor- 



relation function in Eq. (45) results in the dilute sol 



vent regime from the response of individual dipoles to 
the induction field. Moreover, from the recurrence rela- 
tion ( |46[ ), one finds that the next contribution of order 
O ((k^oJ'') is characterized by a total range of 2a, 



C2{Z) 



1 

480 



{ (-z^ - lOz" + 40z^ - 402^ + 12) 61(1 
+ {2-zf9{z~l)e{2-z)Y 



(48) 



with the rescaled distance z — z/a. To conclude, each 
correction term of this expansion corresponding to a 
higher order contribution in the solvent concentration ex- 
tends the range of non-local polarization effects by one 
molecular size. 



Substituting now the relation (45) into Eq. (39), one 



gets for the electrostatic field the following expansion in 
powers of the solvent density, 

E{z) = D{z)-^^iBY,{n,af''Pn{z). (49) 

n>l 
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FIG. 3: (Color online) Dielectric susceptibility function ob- 
tained from Eq. (44 1 normalized by the susceptibility of the 
dielectric continuum approach Xb = psbPo/i for po ~ Qa. — 1 
A kept constant, and different values of a. 



where the lowest contribution to the dielectric response 
is given by the dimensionless polarization field Pi{z) — 
J dz' Ci{z— z')D{z') / {4:TT£Ba) , and the higher order terms 
with n > 1 are obtained from the relation 



Pn{z) 



dz' 



Ci{z-z')Pr,-l{z'). 



(50) 



Considering our preceding discussion, the expansion in 
Eq. ( 49 1 indicates that at the lowest order in the solvent 
density O ((ksi)^), the dielectric screening is induced by 
the response of individual dipoles to the induction field . 

Furthermore, from the recurrence relation dSOl), it fol- 
lows that the next correction of order 



Eq. ( 49 ) with an alternating sign is induced in turn by 



the response of all individual dipoles to the polarization 
field of order O ((/tso)^) , and this correction to the polar- 
ization field positively adds to the induction field. Thus, 
collective effects come into play at the order O ((kso)"'). 
fndeed, the expansion (49) shows that the same self con- 



sistent relation between the response of individual dipoles 
and the resulting modification of the polarization field 
can be extrapolated to higher orders n. We also note the 
apparition of this cooperative behavior at the MF level 
of approximation is rather remarkable. 

We finally show in Fig. [3] that while keeping the dipole 
moment po constant and decreasing the dipole size a, the 
susceptibility function becomes gradually more localized. 
Taking the dielectric continuum limit a — >■ 0, the non- 



locality disappears, and one gets from Eq. (44 1 



lini - z') = XhSiz - z'), 

a— J-0 



(51) 



with Xb = PoPsb/3 the dielectric susceptibility of the PB 
formalism in the region k,;z <^ 1. Thus, in the point 
dipole limit, the dielectric response of the polar medium 
to an external field becomes local, and one recovers the 
dielectric contimmm result E(z) = e^D[z) for the net 
external field. 



p = 0.01 M 




FIG. 4: (Color online) Solid curves : excess solvent parti- 
tion function Eq. ( 52 1 against the rescaled distance from the 



charged interface at vanishing salt density and several bulk 
concentrations. The surface charge is Us — 0.05 e nm~^. Tri- 



angles denote the one-particle density profile of Eq. (All I 
circles are from Eq. (55 1 at the order 0{(^K,sa}^ 
mark the high concentration limit in Eq. (156 



and squares 



B. Inter facial solvent densities 

For low surface charges corresponding to weak poten- 
tials (j){z) < 1, the excess number density that we define 
as 5ksp{z) — ksp{z) ~ 1 follows from Eq. ( 11 ) in the form 



Sksp{z) 



2 



dk cos{kz) 
^Kf + k'^e{k) 



sin(fca) 



ka 



, (52) 



where the function e(fc) introduced in Eq. (21) is given 
by 



(53) 



with the dipolar screening parameter = STrlsPsbQ^ 
and charge structure factor F{x) = 1 — sin(x)/a;. We 
note that for the dipolar geometry and at the physiolog- 
ical solvent concentration psb = 55 M, the bulk dielectric 
permittivity of Eq. (22 1 is given by the Debye-Langevin 
equation Sw = 1 + ^-KiBQ^a^ Psbl'^ — 76.75. 



The excess solvent partition function in Eq. ( 52 ) is dis- 
played in Fig. [4] for various dipolar bulk concentrations. 
First of all, it is seen that at all concentrations p^f,, there 
exists a net solvent excess in the neighborhood of the in- 
terface, and the dipolar attraction is weakened with an 
increase of psb- Then, we notice that beyond the dilute 
solvent regime p^b ?i 1 M, the increase of the bulk solvent 
density also results in a reduction of the range of the in- 
terfacial dipolar attraction. In order to elucidate these 
points, we will consider the opposite limits of dilute and 
concentrated solvents where close form expressions for 
solvent densities can be derived. 
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In the dilute solvent regime KgO < 1, by expanding in 
the region KiZ <^ 1 the integrand of Eq. ( 52 ) in powers of 



K^a, or comparing Eqs. (41| and (45 



between the charge and excess num 
tions 



with the equality 
jer partition func- 



ksciz) = 2Q6k,p{z), (54) 
one obtains for the excess solvent partition function 

5K^{z) = - V(-l)"-^ (^sfl)'""' C,,{z), (55) 



where the functions C„{z) are introduced in Eq. (46) 



First of all, we recognize in the first term of the se ries 
the one particle dipolar partition function Eq. (All I de- 
rived in Appendix |A] This contribution was shown in the 
previous part to result in the polarization field of indi- 
vidual dipoles (the first term of the series in Eq. ([49|). 
This limiting law is reported in Fig. |4]for psi, — 0.01 M. 
Thus in the dilute solvent regime, the interfacial dipolar 
excess extends exactly over the distance a, and the dipole 
behaves as an overall neutral molecule for z > a. 

Then, noting that the function C2{z) of Eq. (481 is 



posit ive for all z, one sees that the next leading term in 
Eq. (551 of order O ((k^c)^) that becomes relevant for 
Kstt ~ 1 (or psb ~ 0.1 M for a = 1 A) corresponds to 
a net reduction of the single dipole partition function. 
This effect is also illustrated in Fig. |4] where we com- 
pare Eq. (55) at the order O ((kjo)^) with the exact MF 
result of Eq. (52 ). The decrease of the local solvent den- 



sity with an increase of the bulk dipole concentration at 
this order results from an intensification of the dielectric 
screening effect that was shown to be induced by the re- 
sponse of individual dipoles to the surface charge. This 
reduction of the polarization density results in turn in 
a decrease of the amplitude of t he polarization field at 
the next order O ((k^o)'') in Eq- (49). These hierarchical 



dipolar response relations can be extrapolated to higher 
orders in solvent concentration by comparing Eq. (55) 
with Eqs. (|49l)-(l50l 



It follows that the modification of 
the single dipole density by collective effects up to the 
order O ((K^a)^") is induced by the response of individ- 
ual dipoles to their own reaction field at to the order 
O ((kso)^""^), and this in turn results in a modification 
of the polarization field at the next order O ((re^a)^"'"'"^) . 
We note that this picture reminds the ionic correlation 
effects in inhomogeneous Coulomb liquids characterized 
by a mutual adjustment of the local ion densities and the 
electrostatic potential with increasing coupling parame- 
ter nais^. 

In the opposite regime of concentrated solvents K^a 3> 



1, the short distance asymptotic limit z/a <^ 1 of Eq. ( 52 ) 



is given by an exponential decay law associated with the 
characteristic decay length kJ^, 



5ksp{z) 



1 
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FIG. 5; (Color online) Fourier transformed dielectric per- 
mittivity function Eq. ( 53 1 for dipoles (solid black curve) , 



quadrupoles (dotted blue curve) and octupoles (dotted red 
curve) for solvent concentration psb = 55 M, and molecular 
size a = 1.0 A for dipoles, a — 0.57 A for quadrupoles, and 
a = 2.68 A for octupoles. The dashed black curve marks the 
dielectric continuum limit e{ka 0) — Sw = 76.75 for the 
dipolar liquid. 



This limiting law, shown in Fig. |4]to perfectly match the 
prediction of Eq. ( 52 ) for psb > 4 M explains the strong 



reduction of the thickness of the interfacial solvent layer 
with an increase of the bulk solvent density. The ex- 
ponential decay indicates that at distances smaller than 
the dipole size, the solvent molecules lead exclusively to 
a charge screening, i.e. they screen the induction field 
D{z) as a concentrated salt solution. 



Effective dielectric permittivity and ion 
densities 



Before considering the variations of the effective di- 
electric permittivity in real space, it is instructive to un- 
derstand the two opposite limits of the Fourier trans- 
formed permittivity function in Eq. ( 53 ) . As displayed 
in Fig. [5] for the biological solvent density pdb = 55 M, 
in the IR limit fca — >■ corresponding to distances much 



larger than the molecular size, the function ( 53 ) tends to 



the bulk permittivity given by the Debye-Langevin rela- 
tion, i.e. e(fc) In the opposite ultraviolet (UV) 
limit fca oo corresponding to the close vicinity of the 
charge source, the permittivity function tends to the per- 
mittivity of the air e(fc) — >■ 1. It is interesting to note 
that the overall shape of this permittivity function resem- 
bles the form of the phenomenological Inkson dielectric 
model |28j. We will investigate below the corresponding 
behavior of the effective permittivity in real space. 



The effective dielectric permittivity profile Eq. (34) is 
displayed in Fig. |6]for various solvent concentrations from 
the dilute to the physiological concentration regime. In 
the dilute solvent regime, by expanding Eq. (34) at the 
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FIG. 6: (Color online) Effective dielectric permittivity for the dipolar liquid at the bulk solvent concentrations (a) = 0.5 M 
(Inset ; psb ~ 0.1 M), and (b) pab = 55.0 M (Inset : rescaled electrostatic field from Eq. (24 1 (blue curve) and Eq. (60 1 (black 



dots)). In the main plots of (a) and (b), solid blue curves are from Eq. (34 1, solid red curves in (a) mark the dilute solvent 



expression in Eq. (571, and the red dashed curve in (b) display the dipolar salt screening regime of Eq. (59 1. The black solid 
curves in (a) (inset) and (b) account for the rotational penalty for dipoles at a rigid interface, (c) Solvent charge density for 
the same model parameters as in (b). (d) Comparison of effective permittivity profiles for dipoles (blue curves), quadrupoles 
(black curves), and octupoles (red curves). In all plots, the solvent concentration is pst = 55.0 M, the solvent molecular size 
a = 1.0 A for dipoles, a — 0.57 A for quadrupoles, and a — 2.68 A for octupoles. 



order 0{{Ksa) ), one obtains for the dielectric permittiv- 
ity profile the close form expression 



£eff{z) ~ 1 



6 



1 



(l-^)%(a-z)|. (57) 



The limiting law (57) is reported in the inset of Fig. [6ja 
(red curve) for psb = 0.1 M. It is seen that the permit- 
tivity increases from the dielectric permittivity of the air 
to the bulk permittivity Sw in a monotonical way over 
one molecular size. This results from the interfacial sol- 
vent charge formation driven by the single dipole-surface 
charge attraction in the air medium. Slightly increasing 
the solvent concentration to psb = 0.5 M (main plot), the 
permittivity curve acquires an oscillatory shape around 
the limiting law (57), and exhibits a peak corresponding 
to a local dielectric increment aX a/2 < z < a. This 
is the density regime where the collective dielectric re- 
sponse mechanism discussed in the previous parts comes 
into play. In Fig.[6jb, it is shown that in a polar liquid at 



the physiological concentration pst — 55 M, the oscilla- 
tory shape of the dielectric permittivity profile becomes 
more pronounced, with the apparition of alternating di- 
electric increment and decrement layers characterized by 
a quasiperiodicity of the order a. Furthermore, at the po- 
sition of the first peak, the local dielectric permittivity 
exceeds the bulk permittivity almost by a factor 2. By 
adding into this picture the rigidity of the interface (see 
Appendix [b] for details), we found that the first dielectric 
increment peak is significantly decreased, and the permit- 
tivity curve is shifted towards larger distances. This re- 
sults from the reduction of the polarization field induced 
by the rotational penalty in the region z < a. 



The oscillations of the background permittivity around 
Ew can be shown to result from the formation of succes- 
sive hydration layers around the charged surface at 2 = 0. 
To this end, we first note that the derivative of the rela- 
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tion ( 43 ) can be written as 



d £off(z) - 1 _ 2psciz) 



dz 



cffiz) 



(58) 



where the solvent charge density is given by Eqs. (52) 



and ( 54 1 . According to the relation ( 58 ) , any reversal in 
the trend of the dielectric permittivity Scs{z) (i.e. any 
minima or maxima) should originate from an alternation 
of the sign of the local solvent charge density. This effect 
is illustrated in the inset of Fig. |6|^c) where we display 
the renormalized solvent charge density against the sep- 
aration distance from the surface. In agreement with the 
number of maxima and minima in the permittivity curve 
of Fig.|6jb), there exists four solvation layers of alternat- 
ing charge over a region of two molecular sizes. 

We illustrate in Fig. [7] the ion number densities 
k±{z) = l^qi(l)o{z) for monovalent ions qi = 1. It is seen 
that the oscillatory shape of the effective dielectric per- 
mittivity in Fig. [6jb results in weak oscillations of the ion 
densities around the PB result k±{z) — 1^ e^^^^ /{qifii). 
Furthermore, the surface dielectric deficiency that in- 
creases the surface potential amplifies the interfacial 
counterion attraction and colon repulsion of the PB the- 
ory. We also show that accounting for the surface rigidity, 
the larger dielectric void in the close neighborhood of the 
interface in Fig. ([6])(b) leads to an amplification of the 
counterion attraction and coion repulsion. However, the 
qualitative behavior of ion densities is not modified by 
the rigidity of the surface. 

At this stage, it should be emphasized that despite 
the simplicity of our linear dipole model and the linear 
MF approximation, the permittivity curve in Fig. is 
able to reproduce qualitatively the shape and the period- 
icity of the transverse dielectric permittivity profiles [29] 
obtained in MD simulations of polar liquids at planar 
interfaces [31 [3]. Furthermore, the overall trend of the 
same dielectric permittivity curve is also in line with the 
result of AFM experiments for water at charged mica 
surfaces, where a rise of the local dielectric permittivity 
profile from e{z) = 4 to with increasing distance from 
the charged surface was observed . This indicates that 
the consideration of the finite size of solvent molecules in 
our microscopic polar liquid model is a key improvement 
over the local theories in order to capture the non-local 
dielectric response behavior of water in real systems. 



It was shown in the previous part IV B that in the 
concentrated solvent regime 3> 1 (i.e. psb 3> 0.1 
M) and at separation distances smaller than the dipole 
size z/a ^ 1, the solvent molecules interact with the 
induction field as a strong salt solution. By computing 
the short distance limit of Eq. (34), one finds that this 



salt screening translates into an exponentially growing 
effective dielectric permittivity. 



effiz) 



(59) 



This limiting law is displayed in Fig. f6|b) by the dashed 
red curve. Substituting now Eq. (|59|) into the rela- 
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FIG. 7: (Color online) Counterion (top) and coion (bottom) 
density profiles of non-polarizable monovalent ions at a pene- 
trable (solid blue curve) and rigid interface (solid black curve) . 
The dashed red curve displays ion densities from the linear 
PB formalism. The model parameters are pib = 0.01 M and 
CTs = 0.05 e nm~^. 



tion ( 33 ) , one obtains for the electrostatic field 



E{z) 



qifj-i 



(60) 



We illustrate in the inset of Figs.[6jc the electrostatic field 
Eq. (pi]) (blue curve) and its asymptotic limit Eq. (60) 



(black dotes) rescaled by the field of the PB formula- 
tion Ep^{z) — e~'^'^ / {q.iPi). First of all, one notices in 
this figure and Eq. ( 60 ) that the PB formalism that can- 



not account for the dielectric screening deficiency on the 
surface underestimates the surface field by a factor e^,. 
Then, we see that the fast drop of the electrostatic field 
to the order of magnitude of Epb{z) is solely driven by 
the dipolar charge screening. This means that at phys- 
iological solvent concentrations, the interfacial decay of 
the surface field is induced by the dipolar salt screening, 
rather than the dielectric screening resulting from the 
preferential orientation of dipoles. However, we note that 
our MF level of approximation neglects image-dipole in- 
teractions, which are expected to weaken the dipolar salt 
screening effect. 



D. Multipolar contributions to non-local dielectric 
response 

We investigate in this part the contribution of the 
multipolar moments of water to the dielectric permit- 
tivity of the liquid. To this end, we incorporated the 
relative multipole/dipole moments of the TIP4P/2005 
water site model into the present theory by computing 
the ratios between the dipolar moment pq = Qa, and 
the quadrupolar and octupolar moments Qq — 
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and = -~Qa^/27 of the solvent molecules depicted 
in Fig. [2ja)-(c). By setting these ratios to the rescaled 
multipole moments Oq/mo = 0.08 A and f^o/A^o = —0.71 
A given in Ref. [24] for the same linear multipole ge- 
ometries as in our model, we obtained the correspond- 
ing molecular sizes. For multipolar solvent molecules. 



the permittivity function ( 53 ) should be evaluated with 



the multipolar charge structure factors and molecular 
sizes given by F{k) = 3 — 8sin{k/2)/k + sm{k)/k and 
a — 0.57 A for quadrupolar molecules, and F(k) — 
10-45 sin(A:/3)/fc-|-9sin(2A:/3)/fe-sin(fc)/fc and a = 2.68 
A for octupoles, where we introduced the adimensional 
wave vector k = ka. 

The shape of the Fourier transformed dielectric per- 
mittivity profiles are displayed in Fig. [5] for dipoles, 
quadrupoles, and octupoles. It is seen that the behavior 
of e{k) is qualitatively similar for quadrupoles and oc- 
tupoles. Namely, the permittivities tend to the air per- 
mittivity for large wave vectors as in the dipolar case, 
which is associated with the dielectric void in the neigh- 
borhood of the charged surface. However, unlike the 
dipolar permittivity curve, both functions exhibit a peak 
corresponding to a region of maximum dielectric screen- 
ing in real space, and converge again towards the air per- 
mittivity for ka — > 0. The latter aspect stems clearly 
from the zero dipolar moment of quadrupolar and oc- 
tupolar molecules. As it will be shown next, this unables 
them to polarize the medium at large distances from the 
charge sources. 

The effective dielectric permittivity profiles associated 
with quadrupolar and octupolar molecules are displayed 
in main plot and the inset of Fig.[6jd). In agreement with 
the wave vector dependence of the Fourier transformed 
permittivity functions in Fig. [5] the dielectric permittiv- 
ities tend to the air permittivity at the charged surface 
and in the bulk limit, with a maximum dielectric screen- 
ing peak in between. This aspect can be explained in an 
intuitive way in terms of the solvent charge density that 
we display in Fig. |6|^c). In this plot, one first notices the 
strong oscillatory behavior of solvent charge densities, 
characterized by sharp kinks located at the separation 
distances between the elementary charges on the solvent 
molecules. Then, for the multipolar liquids, one notices 
the large amplitude of the second negative solvation shell 
following the first positive one at the interface. Hence, 
unlike the dipolar solvent molecules that result in a net 
positive accumulated charge in the interfacial area, the 
first positive charge layer of the quadrupolar and octupo- 
lar liquids are almost exactly canceled by the next neg- 
ative solvation layer. According to Eq 



(43) 



this results 

in a vanishing multipolar contribution to the bulk effec- 
tive permittivity of the medium. Finally, we show in the 
main plot of Fig. that even in the interfacial region, 
the background dielectric permittivity induced by dipoles 
largely dominate the multipolar one. This observation is 
in line with recent MD simulations where the multipolar 
moments of water molecules were shown to weakly affect 
the transverse permittivity of the polar liquid [9j. 



E. Ionic polarizability 

This part is devoted to the effect of ionic polarizabil- 
ity on ionic partitions and the dielectric propreties of a 
dipolar liquid. The charge composition of polarizable 
ions of two species with an equal bulk concentration pi,i 
and electronic cloud radius bp is illustrated in Fig. [2Fd) . 
The elementary charges on the molecules have valency 
e± — =f1 and c± = ±2, with the subscripts ± denoting 
the overall positive and negative molecules in Fig. ^d) . 

For this ionic charge geometry, the susceptibility func- 
tion introduced in Eq. (27) takes the form 



PlPsb 

2a 



{a~\z\) 

2 ^ 



(61) 



+8p.,6,^^exp(^-|^ 



M 

2br, 



Erfc 



M 

26„ 



We note that the function in the bracket on the rhs of 
Eq. (61) behaves as ~ z~^e~^ for z — z/bp ^ 1. This 



indicates that the induced ionic polarizability extends 
the range of the liquid polarizability beyond the solvent 
molecular size a, with a fast decay that obeys a gaussian 
law characterized by the decay length bp. Furthermore, 
for the same charge geometry, the Fourier transformed 
dielectric permittivity in Eq. (21) takes the form 



-A' 

fc2 



sin(A:a) 
ka 



32TT£BPib 



-bik' 



(62) 



The bulk dielectric permittivity that follows from the IR 
limit of Eq. (62) is given by = 1 -I- AniBPoPsb/^ + 
32Tr£BbpPib- Identifying the ionic polarizability as a = 
46p, one finds that the IR limit of Eq. ( 62 ) yields the cor- 
rection from the induced polarizability to the bulk per- 
mittivity derived in Ref. j30| . 



The behavior of the function ( 62 1 is illustrated in the 
inset of Fig. [sjja) at the solvent concentration psb = 55.0 
M, and two values of the bulk ion concentration. In or- 
der to clearly illustrate the contribution from the induced 
polarizability, a considerably large value &p = 5 A was 
chosen. It is seen that the ionic polarizability affects the 
permittivity function mainly at small wavelengths. More 
precisely, for finite polarizability with bp > a, the per- 
mittivity function exhibits a bimodal decay at the wave- 
lengths corresponding to the average fluctuations of the 
electronic cloud radius fci ~ and the solvent size 
^2 ^ a^^. The corresponding behavior of the effective 
permittivity in real space is also shown in the main plot 
of Fig. |8][a) . In agreement with the trend of the function 
e(fc) in Fourier space, the local dielectric permittivity in- 
creases with the induced polarizability. 

In the weak surface charge regime and for the charge 
composition depicted in Fig. [2j^d), the number den- 
sity partition function of polarizable ions follows from 
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FIG. 8: (Color online) (a) Effect of induced ion polarizability on the dielectric permittivity in real space (main plot) and Fourier 
space (inset), (b) Counterion (top) and colon (bottom) densities of polarizable ions for the model parameters pij, — 0.01 M, 
and as — 0.05 e nm~^. 



Eq. ( 10 ) as 



SUMMARY AND CONCLUSIONS 



dk±p{z) = Myjas 



dk cos(fcz) 



(63) 

This density profile is displayed in Fig. [8|Jb) for different 
values of the induced polarizability bp. One sees that in 
the presence of a finite polarizability, the trend of the 
densities in the interfacial area are completely reversed. 
Namely, for weak polariz abilities bp < a, the colon den- 
sity reaches a minimum, and starts to increase towards 
the surface with decreasing separation distance from the 
interface, while counterion density exhibits a concentra- 
tion peak at a characteristic distance, and decreases to- 
wards the interface. For ionic polarizabilities close to the 
solvent molecular size bp ~ a, the interfacial reduction 
of the counterion attraction and colon depletion becomes 
monotonous. 

The surface propensity of colons and exclusion of coun- 
terions with finite polarizability in Fig. [sjjb) clearly re- 
sults from their discrete charge structure. Indeed, the 
ionic polarizability favors the interaction of the negative 
(positive) charge e, on the counterion (colon) with the 
surface charge. This result in a reversal of the ion par- 
tition trends at the interface. Hence, unlike the effective 
permittivity of the liquid in Fig.jsjja), the ion densities are 
substantially affected by the induced polarizability. This 
result disagrees with the conclusion of Ref. [S] where 
the consideration of the ionic polarizability in the point 
dipole limit was shown to weakly affect the ion densities 
in the weak electrostatic coupling regime. This shows 
that our proper treatment of ionic polarizability by ex- 
plicitly considering the extended charge structure of ions 
is crucial. 



In conclusion, we have presented a microscopic theory 
of non-local electrostatic interactions in polar liquids. It 
was shown that unlike previous approaches treating the 
solvent molecules as point dipoles [^ [TBHTB] , our formula- 
tion accounting for the finite size of solvent molecules can 
qualitatively capture the non-local dielectric response of 
polar liquids at charged interfaces. 

In the first part of the article, we derived the field 
theory of the polar liquid composed of linear multipoles 
of finite size, and containing polarizable ions modeled 
as Drude oscillators. From the saddle point solution of 
the partition function, we obtained a non-local Poisson- 
Boltzmann (NLPB) equation. In the rest of the article, 
we investigated the non-local electrostatic interactions 
embodied in this equation within the linear dielectric re- 
sponse regime for polar liquids and ions in contact with 
a weakly charged planar interface. 

In the second part of the article, we introduced a map- 
ping from the microscopic model to the macroscopic for- 
mulation of non-local electrostatics. A key result of this 



part is the expression (43) for the background dielectric 



permittivity of the medium in terms of the accumulated 
polarization charge between the charged interface and 
the liquid. In agreement with MD simulations [31 [HI [TU] 
and AFM experiments of water at charged surfaces [5], 
this relation predicts an interfacial layer associated with 
a reduced dielectric permittivity, and resulting from the 
reduction of the polarization field towards the interface. 

Then, in the third part, we thoroughly analyzed the 
dipolar correlations in the solvent model. We found that 
the non-local dielectric response of the liquid to charge 
sources is driven by a cooperative mechanism resulting 
from the response of the solvent molecules to their own 
polarization field. We also showed that our model can 
qualitatively reproduce the shape and the periodicity of 



14 



the transverse dielectric permittivity profiles obtained in 
MD simulations of water at charged interfaces [3l [9] . The 
fluctuations of the dielectric permittivity around the bulk 
one was shown to result from the formation of successive 
hydration layers of alternating net charge in the interfa- 
cial region. At the next step, we evaluated the contribu- 
tion of the multipolar moments of water to the dielectric 
permittivity of the medium, and found that the multi- 
polar contributions are largely dominated by the dipolar 
one. This observation is in line with the evaluation of 
multipolar contributions to the transverse dielectric per- 
mittivity in MD simulations [31 [S] . 



Finally, we investigated the effect of the induced ion 
polarizability on the ionic partitions at a charged inter- 
face. It was shown that in the presence of an arbitrary 
surface charge, the polarizability completely reverses the 
interfacial ion density predictions of the PB approach, 
resulting in a surface propensity of colons and depletion 
of counterions. This result disagrees with previous works 
based on the point dipole approximation, where a pertur- 
bative correction from the induced polarizability to ion 
densities was observed [3T]. This indicates that the con- 
sideration of the extended charge structure of polarizable 
molecules is crucial. 



Being a first microscopic theory of non-local elec- 
trostatic interactions, the NLPB approach possess 
limitations. First of all, the present theory neglects ex- 
cluded volume effects associated with solvent molecules 
and ions. This complication could be incorporated into 
the theory by imposing a steric Fermi distribution to 
particles [32] or modeling the hard-core repulsions be- 
tween them with a Yukawa potential as in Refs. [33Vl35j . 
Furthermore, we focused in the present work exclusively 
on the linear dielectric response regime of liquids in 
contact with a weakly charged single interface. We wish 
to consider the non-linear effects embodied in Eq. ( 12 ) 



for polar liquids in confined geometries in an upcoming 
article. Then, our discussion on the non-local dielectric 
response of the liquid was based on the MF formulation. 
However, it should be noted that through the field 
theoretic formulation of the multipolar liquid model in 
Eq. ([9|, the present work sets non-local electrostatics on 
a solid theoretical framework. Accompanied with MC 
simulations of the solvent model Eq. Q at interfaces, 
this consistent framework will allow to consider in the 
future non-local electrostatic correlations in inhomo- 
geneous polar liquids in a systematic way. Finally, 
the linear multipole model can be easily generalized 
to different water site models used in MD simulations. 
These extensions will allow direct comparisons of the 
theory with MD simulation results of biological and 
interfacial systems with explicit water. 
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Appendix A: Low density expansion 

This appendix is devoted to the evaluation of single 
dipole densities in contact with a plane located at z = 
and carrying a surface charge Ug. In the dilute liquid 
regime, the expansion of the adimensional grand poten- 
tial JIg = — In Zq in powers of the particle fugacities 
yields the grand potential in the form 



(Al) 



where the statistical average is taken with the gaussian 
Hamiltonian 



dr 



[V0(r)]^ 



ia{r)(j){r) 



(A2) 



Furthermore, the gaussian part of the grand potential 
in Eq. (All is given by f2o = — IhZq, and the reference 



partition function reads 
Zo = j Vcf) e--f^»['^l 
= \/ det(wc) exp 



(A3) 



drdr' 



CT(r)wc(r - r')cr(r') 



Evaluating the field theoretic averages in Eq. (Al), one 
obtains the grand potential in the form 



(r)-V'.(r) 



(A4) 



-A,,e' 



47r 



-Wsi(r)-W's2(r,r+a)-i/'a(r,a) 



where we introduced respectively the ionic and dipolar 
potential of mean forces (PMFs) 

= I dr'dr'V(r')«c(r' - r")g,(5(r" - r) (A5) 

V's(r,a) = j dY'dv"a{r')vc{v' -v")Q[5{r" (A6) 

-<5(r" - r - a)] . 



We recognize in Eq. (A6) the coupling potential in the 



air medium between a single dipole and the fixed surface 
charge. Evaluating the integrals in Eq. (A6| with the 
Coulomb potential 



Vc{y - r') = AttIi 



d^k e^'' ('-'"') 
(27r)3 P 



(A7) 
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one obtains the dipolar PMF in the form 

\z\ \z + a^\ 



(A8) 



We notice that Eq. (A8| is simply the ID interaction 



potential of a finite size dipole in the air medium with a 
constant electric field = l/{Qfis) = 2Tr£B(Js- 
We now note that the dipole density is defined as 



6n 



ip 



Psb 



-Wsi {r)-Ws2 (r,r+a)- V's (r,a) 



47r 



(A9) 



where we accounted for the relation between the fugac- 

ity and the bulk dipole concentration psb — A^e'^ ~. 
Setting the dipolar wall potential to zero and expanding 
Eq. ( A9 1 at the linear order in ipg , one gets the dipolar 



partition function as 

ksp{z) = 1 



'2a 



(AlO) 



Evaluating the integral in Eq. (AlO) with the PMF (A8|, 



one finally obtains the excess dipolar partition function 
in the form 



(All) 



Appendix B: Evaluation of the electrostatic 
potential at rigid interfaces 

We present in this appendix the evaluation of the elec- 
trostatic potential for the polar liquid symmetrically par- 
titioned around a charged rigid interface. The symmetric 
ion distribution around the interface leads to the van- 
ishing ionic wall potential Wi{z) = 0, whereas the sur- 
face rigidity that restricts the dipolar rotations in the 
region \z\ < a can be taken into account by introduc- 
ing the dipolar potential Ws{r, ft) = Wd{a, a + a^), with 
Wd{z, z + Qz) = if z{z+az) > 0, and Wd{z, z+az)) = oo 
for z{z + az) < 0. In the weak potential approximation, 



the MF equation ( 12 ) accounting for this dipolar wall 
potential reads 



in Fourier space. We will thus solve this equation by 
using a perturbative inversion method. To this aim, we 
reexpress Eq. (Bl) in the form 



/a 
[cb{z + az) ~ (f>(z)] 
-a 2a 

^ ~4TTiB [asiz) - XdPsbSksciz)] , (B3) 
where the excess dipolar charge density reads 

Sk^z) = 2Q^e{z)e{a-z) f ^' ^ [d^iz + az) ^ Hz)] 

2Q'e{~z)0{a - \z\) r ^ + az) - cf,(z)] . 
J\z\ 2a 

(B4) 



We note that in Eq. (B3), we introduced the expansion 



parameter in order to keep track of the perturbative 
order. With the use of the electrostatic kernel Eq. ( 16 ) 
that reads for the dipolar charge distribution 



G-\z,z') 



B 



S{z- 



AttI_ 

fPsb J ^ {S{z' - z~az) 



(B5) 



+6{z' ^ z + az) - 2S{z' ^ z)} , 



one can invert the relation (|B3j) and express the potential 
in the form 



(z) = (jjaiz) 



(z), 



(B6) 



where (t>oiz) corresponds to the electrostatic potential of 
Eq. (p3|) for a permeable surface, and the excess potential 



accounting for the corrections from the rotational penalty 
reads 



S^z) 



B J- 



dz'Go{z,z')6ksc{z') 



(B7) 



Furthermore, by expanding the excess potential in powers 
of the parameter Ad, 



(B8) 



substituting this expansion into Eq. (B7| with Eqs. (B4) 



A<j){z) - e^K^,(t){z) + AttIb [us{z) + Psbksc{z)] = 0, (Bl) (|B6j), and identifying the equal powers of A^, one gets the 



with the solvent charge density 



following recurrence relation between the components of 
the excess potential. 



fc,,(z) = 2Q^9iz) 

+2Q^0{-z 



daz 
2a 



[0(z + az) - (f>{z)] 



— min{a.z) 

min(a,|^|) 



— a 



2a 



[(f>{z + az) - (f)iz)] 



(t)n{z) = -{Ksa)^ [ — r(z, z') 

Jo a 

" da. 



(B9) 



X ' [(t>n-liz' - az) - (pn-liz')] , 

la 



(B2) where we introduced the function 



Because the integral boundaries in Eq. (B2) depend on 



the distance from the interface, we cannot solve Eq. (Bl ) 



nz,z')^'rd^^'^^'^^^ (BIO) 
TT Jo [ma) £,„ + q'^e{q) 
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with the rescaled wave vector q = k/a. We note that de- 



riving the expression (B9), we made use of the reflection 
symmetry of the potential 0(— 2) ~ <j){z). 

The steric corrections to the electrostatic potential 
4>n{z) associated with the interface rigidity is evaluated 



numerically from the recurrence relation in Eq. ( |B9[ ) for 
increasing n until numerical convergence is achieved. The 
dielectric permittivity profile in Fig. [6]Jb) (solid black 
curve) was obtained by injecting the converged result into 

Eq. (laal). 
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